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ABSTRACT 

We perform collisionlcss iV-body simulations of 1:1 galaxy mergers, using models which 
include a galaxy halo, disc and bulge, focusing on the behaviour of the halo component. 
The galaxy models are constructed without recourse to a Maxwellian approximation. 
We investigate the effect of varying the galaxies' orientation, their mutual orbit, and 
the initial velocity anisotropy or cusp strength of the haloes upon the remnant halo 
density profiles and shape, as well as on the kinematics. We observe that the halo den- 
sity profile (determined as a spherical average, an approximation we find appropriate) 
is exceptionally robust in mergers, and that the velocity anisotropy of our remnant 
haloes is nearly independent of the orbits or initial anisotropy of the haloes. The rem- 
nants follow the halo a nisotropy - local density slope (j3 — 7) relation suggested by 
lHansen fc Moore] (|2006f ) in the inner parts of the halo, but /3 is systematically lower 
than this relation predicts in the outer parts. Remnant halo axis ratios are strongly 
dependent on the initial parameters of the haloes and on their orbits. We also find 
that the remnant haloes are significantly less spherical than those described in studies 
of simulations which include gas cooling. 

Key words: Methods: iV-body simulations - Galaxies: kinematics and dynamics - 
Galaxies: interactions - Galaxies: haloes 



1 INTRODUCTION 

Simulations of structure formation in the ACDM cosmology 
paradigm consistently produce triaxial dark matter haloes 
with density profiles with p — *■ 00 as r — » 0. This is com- 
monly referred to as a density "cusp". In particular, it is 
found that the spherically-averaged density profile in the in- 
ner halo follows a power law, p oc r m with 70 ~ 1, while in 
the outer parts of t he halo p oc r~ 3 (e .g. iNavarro et al.|[T996l : 
iMoore et aflll998l ; IPwver et all 120031 ). This applies univer- 
sally over all mass ranges investigated. 

The asphericity of dark matter h aloes is also a generic 
prediction of ACDM simu lations (e.g. iDubinski fc Carlberd 
ll99ll ; lAllgood" et al ] |200fJ) , with typical minor-to-major axes 
ratios c/a ~ 0.6 — 0.7 for a Milky Way sized galaxy. However, 
observations of the tidal stream of the Sagittarius dwarf 
spheroidal have been used to argue that the halo of the 
Milky Way is nearly spherical, with c/a < 0.7 ruled out to a 
high degree of confidence for Galactocentric distances in the 
range 16kpc < r < 60kpc (|lbata et allbOOll ^Maicwsk i et ail 
120031 ). This is not necessarily a surprising result, since sim- 
ulations which take into account the effect of gas cooling on 
dark matter haloes find that it results in haloes that are sig- 
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nificantly mo re spherical than those found in pure CDM sim - 
ulat ions (e.g. Katz fc Gunn 199ll : iKazantzidis et alj|20t)3 ). 

lHansen fc MooreT i 2006 ) found that there is also an 
apparently universal relationship between the local den- 
sity slope and velocity anisotropy of simulated dark matter 
haloes. They argued that this could be understood through a 
recognition that the density slope in the tangential direction 
was zero (for a spherical halo), and that the shape of a ve- 
locity distribution is directly dependent on th e density slope 
(whic h is known for simple structures; e.g. lHansen et all 
2005). The anisotropy of the halo is parameterised by 
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and the local density slope as 

„ dlogp(r) 
n ' dlog(r) ' 

lHansen fc Moore] suggested that the anisotropy relates to 
the density slope as 



(2) 



0(7) = 1 -£(1-7/6). 



(3) 



They used data from collapse, merger and cosmological sim- 
ulations to support their hypothesis, and found that the 
best fit to their data come from taking the free parame- 
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ter £ ~ 1.15. lDehnen fc McLaughlin! (|2005l ) showed that as- 
suming this relationship, spherical symmetry, and the re- 
lationship plot oc r~ a where a is some constant (e.g. 
iTavlor fc Navarrdl200ll ). defines a single analytical dynam- 
ical equilibrium for a dark matter halo, independent of a 
with a density profile which closely resembles that found in 

simulations. 

since the work of iToomre fc Toomrd l|l972h on the 
tidal origin of bridges and tails, it has been recog- 
nised that elliptical galaxies could be formed by the 
merger of two disc galaxies. Various authors have 
conducted numerical (JV-body) simulations of collision- 
less mergers, some focusing on the st el lar disc com- 
ponent (e.g. iNeeroponte fc White! Il983l; IBarnesl 1 19881 ; 



iHernquist fc Ostrikerl 1 19921 ; iNaab fc Burkertl 120031), some 
on the dark matter halo (ag. IWhitd Il978l: IVillumse: ' 
19831; iBovlan-Kolchin fc Mai |2004 IKazantzidis et all |200 



Aceves fc Velazquezl l2006; Ka zantzidis e t al. 2006). 
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The studies of White, Villumsen and Boylan-Kolchin & 
Ma, and most other studies of the halo component used ide- 
alised, spherically symmetric systems to represent the galax- 
ies. There are relatively few studies of the remnant haloes 
of mergers between ha l oes co ntaining disc galaxies. 

IKazantzidis et all |2004) looked at the shapes of dark 
matter haloes, primarily considering the effect of gas cool- 
ing, but also reporting simulations of collisionless mergers of 
galaxies with disc components. They stated that the shape 
of the remnant halo depends sensitively on the relative in- 
clination of the discs. lAceves fc Velazquez! (|2006T ) conducted 
simulations of mergers of galaxies with mass-ratios of 1:1, 1:3 
and 1:10, determining that in each case the initial "cu spy" 
density profile was preserved. IKazantzidis et all (|2006l ), as 
part of a larger study of dark matter halo mergers, included 
simulations of binary (1:1 mass-ratio) mergers of haloes with 
disc components. They too found that the presence of a disc 
did not affect the conclusion that "dissipationless mergers 
result in remnants that are practically scaled versions of 
their progenitors". 

All three of these p apers create thei r initial conditions 
using the prescription of lHernouistl l|l993h , which makes the 
approximation of a locally Maxwellian velocity distribution, 
with mean velocity and the dispersion tensor found from 
the moment equations. This method is far from rigorous as 
the true equilibrium velocity distribution can be strongly 
non -Maxwellian. 

INovak et al.l (|2006l ) looked at the shapes of both the 
stellar and dark matter components of the remnants from 
SPH simulations of major mergers. They found that the stel- 
lar remnants were generally oblate, while the haloes were 
prolate or triaxial. The initial haloes in their simulations 
were spherical, and the majority of remnant haloes were 
still relatively close to spherical, having c/a > 0.75, with no 
haloes from any of the 58 simulations having c/a < 0.6. 

The most useful observational work on individual dark 
matte r halo shape s for elliptical galaxie s, iBuote fc Canizaresl 
1 19961 . 1 19981 ) and IBuote et all (|2002l V uses measurements 
of the flattening of X-ray isophotes to place constraints 
on the ellipticity of the dark matter haloes of elliptical 
galaxies NGC1332, NGC3923 and NGC720, under the as- 
sumption that gas rotation has a negligible effect. They 
find that all the haloes are substantially flattened, with 
0.28 < c/a < 0.65, and place tight constraints on the el- 



lipticity of NGC720, which has c/a = 0.38 ± 0.05 for any of 
their oblate density models, and c/a = 0.37 ± 0.04 for any 
of their prolate models. 

No previous simulation studies give a thorough global 
overview of the properties of the haloes of merger remnants, 
covering their structural and kinematic aspects. We perform 
collisionless A-body simulations of 1:1 mass ratio mergers 
of galaxy models consisting of bulge, disc and halo compo- 
nents. A new method for creating a self-consistent equilib- 
rium model of these components is introduced, which avoids 
the use of the local Maxwellian approximation. We examine 
the effects of the initial galaxy orientation, pericentre sepa- 
ration of the mutual parabolic orbit (which we also refer to 
as the "impact parameter" ) , initial anisotropy of the galaxy 
halo, and cusp strength on the density profile, shape, and 
kinematics of the remnant halo. 

In Sections 12.11 and 12.21 we describe the galaxy models 
used in our simulations; in Section 12.31 the various suites 
of simulations we perform are detailed; the density profiles, 
kinematics, and shapes of the remnant haloes are examined 
in Sections 13.11 13.21 and 13.31 respectively. Finally we discuss 
our results and draw conclusions. 



2 GALAXY MODELS AND SIMULATIONS 
2.1 Initial Conditions 

The problem of building an equilibrium disc-bulge-halo 
system is a long standing issue, ad dressed in many 
different ways in the li terature (e.g. iHernquistl 1 19931 ; 
IWidrow fc Dubinski l2005h . Our approach is new, is de- 
scribed in Mc Millan! (120061) , and will be presented in a fu- 



ture paper ( McMillan & Dchncn 2007) . This approach bears 
some relation to that of IBarnesl 1 19881 ). in that the non- 
spherically symmetric component (the disc), is grown adia- 
batically within the spheroid. The process has three stages: 

(i) Creating an equilibrium initial A-body representation 
of the spherically symmetric components of the galaxy (the 
halo and bulge) in the presence of an applied potential field 
corresponding to the monopole (spherical average) of the 
desired disc potential. The distribution function of the halo 
and bulge is that of ICuddefordl l|l99ll ). This distribution 
function has 



/3(r) 



r 2 + I3 rl 
r 2 + r\ 



(4) 



This allows for constant [3 (as r a — > oo) o r for Osipkov- 
Merritt-mode l-like anisotropy (with (3o = 0) (|Osipkovlll979l ; 
lMerrittlll985l ). 

(ii) Evolving the A^-body system whilst growing the non- 
monopole components of the disc potential adiabatically in 
simulation, allowing halo and bulge particles time to relax 
into the full potential of the disc. 

(iii) Populating the disc component with particles. We 
use an im plementa t ion of the disc distribution function /new 
defined in iDehnenl l|l999l ). see Equation [7] of this paper. 

For the halo we use a spherically-symmetric truncated 
generalised Navarro et al. (1996, NFW) like profile: 
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planer component of the disc distribution function as 
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Figure 1. Rotation curve in the plane of the disc for our stan- 
dard 70 = 1 galaxy model, shown for both the initial conditions 
(upper), and after being allowed to evolve in isolation for 300 time 
units (approximately the time it takes the galaxies to reach peri- 
centre in our merger simulations — lower) The solid line is the net 
rotation curve, also shown are the decomposed contributions from 
the disc (dotted), bulge (short- dashed), and halo (long- dashed). 
The evolution seen is dominated by that in the inner ~ 2R d of 
the disc component. This is caused by a bar instability in the disc 
which redistributes its mass. This, in turn, has some effect on the 
halo at its innermost radii (see Section l2.3t . 



where p c is a scale density, is the halo scale radius, 
r t is the halo truncation radius, and 70 describes the in- 
ner slope of the density profile. As r^O, the halo density 
Ph ocr -70 . We truncate at a radius significantly larger than 
the halo scale radius, so in the outer halo the density goes 
as ph oc r _3 sech(r/r t ). Some truncation is necessary as the 
untruncated profile (rt —* 00) is infinitely massive. 

The disc is defined as having a surface density profile 
which is exponential in the cylindrical radius R, and a verti- 
cal (z-component) structure modelled by isothermal sheets 



p d (R, z) = 



M d 



-exp 



-AW (JL 
Rd ) \ zd 
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where Md is the total disc mass, R d is the disc scale radius, 
and Zd the scale height. 

To find the distribution function of the disc, we as- 
sume that it can be decomposed into its components 
in the plane of the disc, and perpendicular to it, i.e. 
fd — fd(E p , L z , E z ) where the "vertical(z-component) en- 
ergy" E z = \yi + &(R, z) — $(-/?, 0), and the "planar en- 
ergy" E P = E-E z =\(vl + v 2 e ) + <S>(R, 0). 

If we define Re p as the radius of a circular orbit in 
the disc midplane with planar energy E p , we can write the 



fd,p{E p , Lz) — 

n(R Ep )X'{RE P ) 

7T k(Re p )ct'£(Re p ) 



exp 



Q(Re p )[Lz— Lz, c (E p 



(7) 



where L Z}C (E) is the angular momentum of a circular or- 
bit in the disc midplane with planar energy E p , £1 is the 
circular frequency of that orbit, and k is the epicycle fre- 
quency. Z'(R) and a'i(R) are defined such that the true 
surface density £(i£) and radial velocity dispersion o%{R) 
of the iV-body representation are those desired, to within an 
appropriate degree of accuracy. 

We model the bul ge as non-rotating , initially spherically 
symmetric and with a lrlernquisj (|l990h density profile. 

M b 
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where Mb is the mass of the bulge, and rb is its scale radius. 

We performed a suite of simulations in an effort to 
investigate the effect of the properties of the progeni- 
tors and of the encounter on the halo of the merger 
of two equal-mass galaxies. The galaxy models were de- 
signed to loosely resemble the Mi lky Way as described by 
iKlypin. Zhao, fc Somerviile] (|2002f l. 



2.2 Numerical Miscellanea 

In all our models we choose units such that the disc scale 
length, Rd — 1, disc mass Aid = 1, and the constant of grav- 
ity G — 1. The disc scale height was chosen to be Zd = 0.1, 
the velocity dispe rsion of the disc was defined such that 
the iToomrd (| 1964T ) stability parameter Q = 1.2 at all radii. 
The bulge mass Mb =0.2, and bulge scale length rb = 0.2. 
Scaling these values to the Milky Way, taking Rd = 3.5kpc, 
M d + M b = 5 x 10 10 M Q gives a time unit ~ 14Myr. We 
chose to have the scale radius of the halo — 6, the trun- 
cation radius r t = 60, and the halo mass Mh = 24. In sim- 
ulations with 70 = 1.0, 79% of this total mass is within 
the truncation radius. The rotation curve for this model is 
shown in Figure [1] 

The stellar components were populated with 150 000 
equal mass particles (i.e. 125 000 in the disc, 25 000 in the 
bulge) with a smoothing length e s teiiar = 0.02. The halo 
component was populated with 750 000 equal mass parti- 
cles with a smoothing length e^ a i = 0.04. That corre- 
sponds to each halo particle being 4 times more massive 
than a stellar particle. A small number of simulations with 
4 times more particles were performed for comparison. These 
demonstrated that the numerical resolution used was suffi- 
cient. 

The iV-body simulations were performed using the pub- 
licly avail able iV-bod y code gyrfalcON, which is based on 
Dehnen's (|2000l : l2002h force solver falcON, a tree code with 
mutual cell-cell interactions and complexity 0(N). The 
equations of motion were integrated using the familiar leap- 
frog integrator with minimum time step 2 -7 and a block- 
step scheme allowing steps up to eight times larger. Individ- 
ual particle time steps were adjusted in an (almost) time- 
symmetric fashion such that on average 
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with "l>i and a,i the gravitational potential and acceleration 
of the ith body. With these parameters, energy was con- 
served to within 0.1% over the full time span (1000 time 
units) in a typical simulation, approximately corresponding 
to a Hubble time. The time span was chosen such that there 
is sufficient time for the remnant to reach a dynamical equi- 
librium. 

In all simulations the equal mass galaxies were placed on 
a mutual orbit in the x-y plane. In the vast majority of cases 
these orbits were parabolic and corresponded to one that 
would have a pericentre separation of 8 Rd if the galaxies 
were point masses; we refer to this distance as the "impact 
parameter" d of the merger. In some cases the galaxies were 
placed on orbits with smaller pericentre separations, or on a 
radial orbit with zero net energy (i.e. a parabolic orbit in the 
limit where impact parameter d — ► 0). The galaxy centres 
were initially separated by 200-Rd- 

2.3 Suites of simulations 

We perform four suites of simulations of equal-mass mergers, 
to examine the effects of varying different parameters on the 
mergers. The suites vary in 

(i) Orientation of the disc galaxy. 

(ii) Pericentre separation of the galaxies' mutual orbit. 

(iii) Initial velocity anisotropy of the halo component. 

(iv) Cusp strength of the haloes (70 in Equation [SJ) . 

In addition, we ran control simulations, which were 
stand-alone simulations of the same galaxy models, run for 
the same length of time, with no merger. The galaxy disc 
is bar-unstable, which in these simulations causes a steep- 
ening in the density profile of the inner parts of the disc, 
which also causes a steepening in the halo profile in the 
inner ~ 0.2rh. This steepening is commonly observed in 
simulations of bar red galaxy evolution (reviewed e.g. by 
I Athanassoula! [2004) . The effect of this evolution on the ro- 
tation curve of the galaxy can be seen in Figure [T] The 
anisotropy of the halo velocity distribution is essentially un- 
changed by the stand-alone evolution. 

2.3.1 Suite 1: Orientation 

Our first suite of simulations consisted of 16 mergers which 
vary only in the orientation of the discs with respect to one 
another, and to the angular mome ntum vector o f their mu- 
tual orbit. Following the example of lBarnesl (|l988i). we define 
four different orientations for the spin vector of each galaxy 
disc, corresponding to pointing towards the four vertices of 
a regular tetrahedron. The orientations of disc galaxies in 
mergers are usually described in terms of the disc inclination 
relative to the orbital plane i, and the argument of pe ricen- 
tre uj (See also Figure l2l or iToomre fc Toomrd \l9m . The 
inclinations (i\ &i2), and arguments of pericentre (uii &W2) 
for the suite of mergers are shown in Table 1. 

These are the same initial orientations that 
were used in t he c ollisionless merger simulations of 
Naab fc BurkertJ (120031) a nd the SPH merger simulations of 
Naab. Jesseit. fc Burkert] |2006T) . These studies detailed the 
significant effects that varying the initial disk orientations 
has upon the properties of the luminous components of the 
merger remnants. 



Cross-section Plan view of orbital plane 

(at pericentre) 

Angular momentum Mine of nodes 

vector of the mutual orbit \ 




Figure 2. Diagram showing the definitions of inclination i, 
and argument of pericentre u> for the mergers. By convention, 
-180° < i ^ 180°, -90° < u < 90°. The "line of nodes" is the 
intersection of the plane of the galaxy disc with the orbital plane. 
The cross section view (left) is that from a point in both the 
plane of the orbit (shown as the horizontal), and the plane of the 
galaxy disc (shown at an angle i to the horizontal). The plan view 
is from above the orbital plane, looking down. 

In all cases the initial halo has an inner density slope 
70 = 1.0 (Equation [5]) , and an isotropic (/3 = 0) velocity 
distribution. The galaxies were put on a parabolic orbit with 
a pericentre separation of 8Rd- 

TABLE 1 



Disc 1 Disc 2 



# 


h 




ii 


UJ-2 


1 








180 





2 








71 


30 


3 








71 


-30 


4 








71 


90 
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-109 


-60 


180 
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-109 


-60 


71 


30 
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-109 


-60 


71 


-30 
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-109 


-60 


71 


90 
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-109 





180 
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-109 
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30 
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-109 
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-30 
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-109 





71 


90 


13 


-109 


60 


180 





14 


-109 


60 


71 


30 


15 


-109 


60 


71 


-30 


16 


-109 


60 


71 


90 



2.3.2 Suite 2: Impact parameter 

We performed two sets of simulations with galaxy models 
identical to those used in Suite 1, varying the separation 
at pericentre (impact parameter) d between 8i?d and (i.e. 
radial impact) in increments of 27?d- The two sets of simula- 
tions differ in the galaxy orientations, which correspond to 
those of simulations 4 & 7 in suite 1 (Table 1) . 

2.3.3 Suite 3: Halo anisotropy 

Our third suite of simulations consisted of mergers of galax- 
ies with varying initial halo velocity dispersions and were 
performed to examine the effect that has on remnant prop- 
erties. The haloes had identical density profiles to those 



The haloes of merger remnants 5 



varying orientation 



varying impact parameter 



Merger remnants 

Stand-alone 

ystand alone density 




Merger remnants 

Stand-alone 

2*Stand-alone density 




Figure 3. Suite 1: Spherically averaged density profile (lower) 
and density slope profile (7, see Equation [2] upper) for remnant 
haloes from the 16 simulations in suite f (solid lines). Also shown, 
for comparison, is the density profile and slope of the halo of the 
equivalent stand-alone simulation (dotted), and the same stand- 
alone density profile with the density doubled - i.e. what one 
would expect if the galaxies were simply stacked on top of one 
another (dashed). 

in suites 1 and 2, but rather than an isotropic veloc- 
ity distribution, they were initialised with constant ra- 
dial anisotropy (/3 = 0.35 throughout the halo), constant 
tangential anisotropy (/3 = —0.4 throughout the halo), or 
Osipkov-Merritt-model-like anisotropy (Equation [4} with 
r a = 24 = 4rh. Simulations were performed with galaxies 
orientated as per orientation 7 (Table 1), and with impact 
parameters of 8, 3, and ORd- 

2.3.4 Suite 4: Cusp strength 

We ran a series of simulations with varying halo cusp 
strengths. The inner slope of the density distribution (70 
in Equation [5]) varied between 70 = 0.1 and 70 = 1.6. In 
all cases we retained the parameters — 6, r t = 60 and 
Mh = 24. The galaxy discs had the same orientations as sim- 
ulation 4 in Table 1. The haloes all had isotropic velocity 
distributions, and the mergers trajectories all had impact 
parameter d — 8Rd ■ 



3 RESULTS 

3.1 Halo density profiles 

The haloes of the merger remnants are not spherically sym- 
metrical (see Section[3]3} , however insight, and a comparison 
with previous work, can be found from looking at spherically 
averaged density profiles, and local density slopes (Equa- 
tion [2j). 

Figure [3] shows the density profiles and slopes for the 



Figure 4. Suite 2: Spherically averaged density profile and den- 
sity slope profile for remnant haloes, as in Figure[3] for simulations 
with impact parameters d = 8, 6, 4, 2, ORd (solid lines). 



16 merger remnants from suite 1 (Section 12.3. lfl . with the 
equivalent for the stand-alone simulations. The disc orienta- 
tions have no effect on the density profile, even in the regions 
where they dominate the galaxies' density. The shape of the 
remnant density profile is very similar to that of the initial 
haloes, with an almost identical scale length. There are, how- 
ever differences in the outer parts: ~ 6% of the halo mass 
becomes unbound during the merger, and the mass within 
10rh is ~ 50% of the total mass of the original haloes (as 
compared to 79% in each halo initially). Hence the density 
within the inner ~ lOr^ of the merged haloes is closer to 
that of one of the initial haloes than to double that (as it 
would be if the merger was equivalent to simply "freezing" 
the haloes, and placing them on top of one another). 

The simulations of suite 2 show that varying the impact 
parameter of mergers also has little effect on the shape of 
the density profile (Figure |4|. Results for suite 4 showed 
the same trends for all cusp strengths. In the interests of 
brevity, only the results for 70 = 0.1 and 70 = 1.4 are shown 
(Figure [5]l . 

Results from suite 3 do show a slight trend in remnant 
cusp strength with initial halo anisotropy. The trend is ob- 
servable for all orbital separations simulated, but only re- 
sults from simulations with orbital separation d = 8Rd are 
shown in Figure [6] for clarity. Haloes which have initially 
constant tangential anisotropy have sharper remnant cusps 
than initially isotropic haloes, while haloes with initial con- 
stant radial anisotropy have weaker remnant cusps. Giving 
the haloes Osipkov-Merritt anisotropy does not alter the 
cusp strength. Compared to the isotropic case, cusp strength 
- measured in terms of l y(r) at small radii - is ~ 0.1 greater 
for the P = 0.35 case than for the isotropic case, and is ~ 0.1 
less in the j3 = —0.4 case. This trend is not seen in the stand- 
alone profiles, which are nearly identical in all three cases. 
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Figure 6. Suite 3: Spherically averaged density profile and 
density slope profile for remnant haloes from d = 8R^ merg- 
ers in which the initial haloes had constant radial anisotropy 
(fl = 0.35: solid line); constant tangential anisotropy (/9 = —0.4: 
short- dashed); Osipkov-Merritt anisotropy (r a = 24: long-dashed) 
Also shown are the density profile and slope for the equivalent 
mergers with initially isotropic haloes, as in all other simulations 
dotted). 



Figure 7. (3 as a function of radius for the simulations of suites 1— 
3 (top to bottom). Error bars shown are typical of all simulations 
in respective plot, found using Equation 1101 In the bottom plot 
remnants of mergers with haloes with different initial anisotropics 
are represented by different line types. Dotted lines plot fi for 
remnants of initial halo /3 = 0.35 mergers, dashed lines for initial 
halo /3 = —0.4, solid lines for the Osipkov-Merritt haloes. 



3.2 Halo velocity distributions 

We investigate the velocity distribution of the haloes by 
looking at (3 as a function of radius, and as a function of 
the local density slope 7()"). (3 is determined for particles in 



spherical shells, with the median particle radius for the shell 
being the value plotted in Figures [7] & [5] 

Error estimates, as shown in Figures[7Jand[S]were found 
for the spherical shell of particles through the standard prop- 
agation of errors formula, such that the error A/3 is given 
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Figure 8. /3 as a function of radius for the simulations of suite 
4. Error bars shown are typical. Simulation results shown are for 
haloes with density cusps of strengths as indicated. 

by 

(A/3) 2 = 

(||)W + (|§)W + (||)W, 

(10) 

The remarkable aspect of the observed 3 profiles of the 
remnant haloes is its near independence from the parameters 
varied in suites 1-3 (Figure [7] top to bottom). In all cases 
examined 8 ~ 0.1 at small r, and increases to/3~0.2 — 0.3 
at r ~ 10rv This shows that while the remnant halo has a 
clear "memory" of the density profiles of the initial haloes, 
it has little memory of the velocity distribution of the initial 
haloes. There is some difference observable between the 3 = 
0.35 and (3 = —0.4 cases, perhaps most noticeable in the 
range 0.5rh < r < r\, but it is vastly smaller than the 
difference in the initial halo conditions. 

One would also naively expect the impact parameter of 
the merger to have a significant effect, but it does not show 
any greater effect on the (3 profile than the orientation of 
the galaxies does. The independence from disc orientation is 
unsurprising at large radii, but also continues to small radii, 
where one might have expected the disc orientation to play 
a role. 

The (3 profiles from the simulations of suite 4 are shown 
in Figure [5] As r — > 0, the value of (3 tends to for the less 
strongly cusped initial conditions (70 = 0.1 — 0.6). As the 
cusp strength increases, the radial anisotropy of the inner 
parts of the halo increases. At larger radii the anisotropy 
tends towards f3 ~ 0.2 — 0.3 for all cusp strengths, with 
the less strongly cusped haloes tending to have higher 
anisotropy. 

Almost all of the (3 profiles share a similar non- 
monotonic shape. (3 increases at small r, reaching a peak 



at r ~ 1 — 1.5rh, then falls off, reaching a local minimum at 
r ~ 2 — 3rh, before rising again. This apparent oscillatory 
variation of [3 is not associated with any tidal arms, shells, 
or any other observed features of the particle distribution; 
it is also not seen to disappear or vary to any great extent 
with time. A similar 3 profi le can be seen in Figure 11 of 
iBovlan-Kolchin fc Mai |2004). Further work is required to 
det ermine its cause and imp ortance. 

lHansen fc Moord 1120061) argued for a universal relation- 
ship between local density slope 7 and velocity anisotropy B, 
given in Equation [3] which we plot in each of the graphs in 
Figure [5] The data fit the suggested relationship for 7 < 2, 
but do not fit for 7 > 2, where 3 is systematically lower 
than the fitting function. While it is the outer parts of the 
halo which have 7 > 2, and the dynamical time is longer here 
than in the inner parts of the halo, this is not simply because 
the outer halo has yet to reach an equilibrium. Simulations 
which ran for twice as long (until t = 2000 = 28Gyr scaled 
to the Milky Way) showed identical results. The 3 = 0.35 
initial conditions produce remnants which are a close fit to 
the suggested relationship, but the 3 = —0.4 initial condi- 
tions produce remnants which have systematically slightly 
lower values of 3 for the same 7. The fitting formula equa- 
tion|3] with the smaller value for the free parameter £ ~ 1.05 
(plotted as a dotted line in figure [9} does appear to provide 
an upper limit for 3. 

We observe that the haloes' net rotation is very small; 
we compare to the velocity dispersion and find that v rot /a < 
0.1 in the inner regions in all cases. There is slightly more 
rotation at large radii, but still v TO t/<? < 0.2 in all cases. 
Bulk rotation is also near negligible, with the axes of the 
bulk distribution remaining near fixed over a long period of 
time in all cases. 

3.3 Halo Shape 

We calculate the overall shape of the remnant halo through 
an iterative procedure based on diagonalizing the moment 
of inertia tensor 

— ^ ITT'afijOtfjjOtj (H) 

where the sum is over all particles considered, and ri >a is 
the i coordinate of the ath particle. We start with a spheri- 
cal window, centred at the densest point of the halo, which 
contains 50% of the mass of the halo, then determine and di- 
agonalize its moment of inertia tensor. This determines the 
principle axes and the axis ratios of an ellipsoidal window 
which is scaled to contain 50% of the halo mass. This is then 
repeated to find new ellipsoidal windows until the results 
converge such that the volume of the window is conserved 
to within a fractional difference of 10 -3 between iterations. 
From the ellipsoidal window we determine the axis lengths 
a, b & c with a b ^ c, and thus the intermediate (b/a) and 
minor (c/a) axis ratios. We can then define an ellipsoidal 
"radius" 

C-f^f^f- (12) 

The question of whether an ellipsoidal body is prolate, oblate 
or triaxial can be quantified through the parameter 

T= (a 2 -b 2 )/(a 2 - c 2 ). (13) 
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Figure 9. Anisotropy parameter /3 plotted as a function of the local density slope 7 for mergers from suite 2 (left), suite 3 (centre), 
and suite 4 (right), fi and 7 are determined for each simulation at many points with the radial range O-lr^ < r < lOr^. The lines drawn 
in each plot correspond to equation(3] with £ = 1.15 (solid line) and £ = 1.05 (dotted). In the central plot, filled squares correspond to 
mergers with j3 = 0.35 initial conditions, open squares are from the = —0.4 initial conditions, and crosses for the Osipkov-Merritt halo 
initial conditions. 
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Figure 10. Minor (c/a) and intermediate (b/a) axis ratios for the 
simulations of all suites of simulations (as labelled). In all figures 
the solid diagonal indicates the position on the plot of a com- 
pletely prolate halo. The dashed lines are of constant triaxiality, as 
indicated. For suites 2 (top-right), 3 (bottom-left) and 4 (bottom- 
right) the different symbols represent different simulation param- 
eters as labelled. The three points for each initial anisotropy (suite 
3) correspond to simulations with impact parameters d = 0,3,8 
with, in each case, b/a^Q < b/a^—3 < b/a^—g. 



3.3.1 Trends in halo shapes 

The graphs of Figure [lOl which show the ratios c/a and b/a 
(determined as described above) for simulations from all of 
our suites, tell us the following things: 

(i) Remnant haloes are significantly non-spherical (c/a < 
0.55 in all cases shown). They are generally prolate, or at 
least tending towards prolate, in our simulations. 

(ii) The triaxiality of the halo at this scale ~ 10>h for 
outermost particles) is barely affected by the initial inclina- 
tion of the disc. 

(iii) More radial mergers (i.e. mergers with lower values 
of d) tend to have more prolate remnant haloes. This is true 
whether or not the initial haloes are isotropic. 

(iv) Tangential anisotropy in the initial halo tends to pro- 
duce haloes that are more prolate, and have a higher value of 
c/a than in the isotropic case. Conversely, radial anisotropy 
in the initial halo tends to produce haloes that are more 
triaxial, and have a lower value of c/a than in the isotropic 
case. 

(v) In the range 0.4 < 70 < 1.6, an increasing cusp 
strength tends to produce remnants with a greater value 
of b/a. This is not simply due to the decreasing size of the 
window containing 50% of the mass in haloes with increasing 
cusp strength, and is also seen if the volume of the window 
is defined such that it is the same for all cusp strengths. 

Like lNovak et all (2006) we find that the minor axis of 
the stellar remnant and the major axis of the remnant halo 
are nearly always close to perpendicular. The major axis of 
the halo is always close to the plane of galaxies' original 
mutual orbit, and the minor axis of the stellar component 
is nearly always close to being perpendicular to that plane. 



We use this metho d to enable a like-fo r-like comparison with 
the simulations of lNovak et al.l l|2006h , though starting with 
a spherical window can cause bias. 



3.3.2 Ellipticity profiles 

We wish to examine the effect of the orientation of the 
disc components upon the ellipticity of the remnant halo 
as a function of radius. In an effort to do this, we find the 
density at each particle by a "nearest neighbours" analysis 
l|Casertano fc Hutlll985l l. divide the particles into shells by 
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spherical shells 




both disks in merger plane 
one disk in merger plane 
neither disk in merger plane 




Figure 11. Minor axis ratios as a function of radius r (upper) 
or ellipsoidal radius £ (lower) for the simulations of suite 1. The 
dotted lines corresponds to the simulation in which galaxy discs 
are both oriented in the merger plane, the dashed lines correspond 
to the simulations in which one of the discs was oriented in the 
merger plane. 



density, and determine the axis ratios and median ellipsoidal 
radius of each shell. The axis ratio of the particles in a shell 
is found from their moment of inertia tensor (Equation 1 
For comparison to iKazantzidis et all (|2004h we also divide 
the merger remnant halo into spherical shells, and find the 
axis ratio of the particles in those shells. We defined the ra- 
dius of the spherical shells as being the median radius of the 
particles in it. The two approaches produced qualitatively 
similar results, though, as I Athanassoulal (|2006l ) shows, using 
a spherical window introduces a bias towards larger values 
for c/a and b/a. In Figure [TT] we plot the results from the 
both analyses. 

Considered in terms of c/a at the innermost point deter- 
mined, the flattest merger remnant comes from the co-planar 
merger, and the next four flattest remnants are from mergers 
with one of the discs in the orbital plane. The difference in 
minor axial ratio is relatively small. Using the density shells 
values, at the innermost point for the co-planar merger rem- 
nant c/a is ~ 0.62, for mergers with one galaxy disc in the 
orbital plane it's ~ 0.67, and for mergers with neither disc 
in the orbital plane it's ~ 0.72. The simulation with orien- 
tation 16 (see Table 1) produces a notably more spherical 
remnant than all the other simulations (c/a ~ 0.83). Be- 
yond r ~ 1.5rh the difference in c/a is negligible. There is 
no trend apparent in the intermediate axis ratio, b/a. These 
trends are exactly as one would expect. 

It is reasonable to ask whether the fact that these haloes 
are so far from spherical invalidates the approach in Sec- 
tion 13.11 of using a spherically averaged density profile. In 
an effort to investigate this, we created density profiles based 
upon the density shells, found through the method described 
above. These profiles, while suffering somewhat from in- 



creased noise, showed the same behaviour as that described 
in Section [3]T] This indicates that the use of spherically av- 
eraged profiles is sufficient and appropriate for the analysis 
of the haloes of merger remnants. 



4 DISCUSSION 

We have performed a number of suites of simulations, cov- 
ering a number of different merger parameters. Across all 
these simulations we have seen that the halo cusp strength 
is extremely robust against changes caused by major merg- 
ers, even when there is a centrally concentrated stellar com- 
ponent. This is in kee ping with previous resu lts from halo 
only simulations (e.g. IKazantzidis et al"1l2006f) ; simulations 
with disc compon ents, but initialised using a Maxwellian ap- 
proximatio n (e.g. lAceves fc Velazquezll2006l ) and analytical 
arguments (IPehnenlbOolf r Thus mergers, at least of gas-less 
equal mass galaxies, can not solve the persistent discrepancy 
between observations of nearby galaxies, which imply that 
galactic dark matter haloes have a density profile with a flat 
core, and the cosmological standard model, which predicts 
that haloes should have a cusp. 

In contrast to the density profile, the velocity anisotropy 
of the halo has very little "memory" of the initial conditions 
of the merger. In all cases examined, the velocity anisotropy 
of the rem nant followed the r elation ship (Equation [3j sug- 
gested by lHansen fc Moorei (^OOd ) relatively closely for 
7 < 2, but [3 was systematically lower than predicted for 
7 > 2. That these results are nearly independent of initial 
anisotropy is very encouraging. The overwhelming major- 
ity of merger simulations to date have been performed with 
isotropic haloes, while the haloes found in cosmological sim- 
ulations are typically radially anisotropic. A strong depen- 
dence on anisotropy would suggest that the approximation 
of isotropy in the halo was a poor one, which would have 
cast doubt upon results from simulations using it. The ve- 
locity anisotropy of the halo is clearly a less robust property 
than its density profile. 

Given t he shape of the dark matter haloes observed by, 
for example iBuote et al.l (|2002h . the results of Section 13.31 
are of particular importance. Axis ratios c/a in the range 
0.37 ± 0.04 (as seen for NGC720) are seen for all initial 
halo anisotropies except (3 — —0.4, in cases with suffi- 
ciently small impact parameters. Si mulations which i ncor- 
porate the effects of gas physics (e.g. lNovak et al1l2006l ) find 
remnant haloes which are significantly more spherical. It 
is likely, therefore, that galaxies with such highly flattened 
haloes were formed from near-radial trajectory, very gas- 
poor mergers in which collisionless dynamics dominate. 

That the remnant halo is closest to completely pro- 
late for mergers with low impact parameters is unsurpris- 
ing, and w as recognised in the literature as far back as 1983 
Villumscn in the case of completely spherical initial models. 
In the d — case, the major axis of the remnant is in the 
same direction as the initial motion of the galaxy centres. 

The dependence of remnant shape on initial anisotropy 
of the haloes is far less trivially understandable, but is likely 
to b e related to the radial-orb it instability (e.g. iBarnesI 
1 19851 ; iDeionghe fc Merrittl I1988T I. Spherical systems with 
large numbers of radial orbits are known to be unstable to 
deformation towards a barred or triaxial shape. The radial- 
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orbit instability is also k nown to reduce the cent ral con- 
centration of models (e.g. iMerritt fc Aguilarl ll985T ), which 
could explain the slight decrease in cusp strength seen in 
the P — 0.35 model remnants (Figure [BJ. 

It is known that a steep cusp in a triaxial model causes 
orbit-s cattering which can ha ve a significant effect on its 
shape (jValluri fc Merrittllf 9981 ). reducing the triaxiality. The 
effect of increasing the cusp strength of our haloes is to make 
the remnants more spherical (less prolate, i.e. increasing 
b/a), but this means that the most strongly cusped haloes 
are more triaxial than the least strongly cusped ones. 

Increasing the cusp strength of the halo puts more of 
the mass of the halo at radii smaller than the impact param- 
eter of the merger. This is similar to the effect of increasing 
the impact parameter of the merger, which we have shown 
increases b/a for the remnant. It is likely that this explains 
the trend in remnant shape with cusp strength. 
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